#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <random>
#include <type_traits>
#include <vector>
#include <stdio.h>
#include <omp.h>
#include "kernels.cuh"
#include <assert.h>
#include <stdlib.h>
#include <math.h>

// Include CUDA runtime and CUFFT
#include <cuda_runtime.h>
// for benchmark
#include <cufft.h>

// Include pocketfft
# include "pocketfft.h"

#pragma once

// Decimation in Time FFT cuda version, used to perform dft of bluestein's algorithm
// Performs bit-reverse-copy before transform
// N must be power of 2
struct DitFFT_cuda {
    static constexpr char Name[] = "DitFFT_cuda";
    const std::size_t N;
    const int L;
    cuFloatComplex* cuda_W;
    Comp* W;
    cuFloatComplex* d_X;
    cuFloatComplex* d_Y;
    size_t* d_I;
    std::size_t* I;

    DitFFT_cuda(std::size_t N) : N(N), L(__builtin_ctzll(N)) {
        using namespace std;
        // initialize W
        auto w = _2Pi / (Real) N;

        // allocate memory 
        cudaMalloc(&cuda_W, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_X, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_Y, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_I, N*sizeof(size_t));

        I = (std::size_t*) malloc(N*sizeof(std::size_t));
        W = (Comp*) malloc(N*sizeof(Comp));

        init_w_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(cuda_W, w, N);

        // initialize I
        
        for (size_t i = 0, j = 0; i < N; ++i) {
            I[i] = j;
            for (size_t k = N >> 1; (j ^= k) < k; k >>= 1);
        }

        cudaMemcpy(d_I, I, N*sizeof(size_t), cudaMemcpyHostToDevice);

        for (size_t i = 0; i < N; ++i)
            W[i] = exp(Comp(0, w * (Real) i));
    }

    ~DitFFT_cuda(){
        free(W);
        free(I);
        cudaFree(d_X);
        cudaFree(d_Y);
        cudaFree(cuda_W);
        cudaFree(d_I);
    }

    void dft(Comp* Y, const Comp* X) const {
        using namespace std;
        cudaMemcpy(d_X, X, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        
        bitrev_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_Y, d_X, d_I, N);

        for(int i = 0; i < L; ++i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditfft_kernel<<<gridDim, blockDim>>>(d_Y, cuda_W, i, N);
        }
        
        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }

    void idft(Comp* Y, const Comp* X) const {

        using namespace std;
        cudaMemcpy(d_X, X, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        
        bitrev_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_Y, d_X, d_I, N);

        double tt = omp_get_wtime();
        for(int i = 0; i < L; ++i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditdifidft_kernel<<<gridDim, blockDim>>>(d_Y, cuda_W, i, N);
            Check_CUDA_Error("267");
        }

        complex_div_real<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_Y, N);
        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);

    }
};

// Forward transform uses DIF while inverse transform uses DIT, used to perform idft of bluestein's algorithm
// For convolution use, no bit-reverse-copy performed, performed in-place
// N must be power of 2
struct DifDitFFT_cuda {
    static constexpr char Name[] = "DifDitFFT_cuda";
    const std::size_t N;
    const int L;
    cuFloatComplex* cuda_W;
    cuFloatComplex* d_Y;
    Comp* W;

    DifDitFFT_cuda(std::size_t N) : N(N), L(__builtin_ctzll(N)) {
        using namespace std;
        auto w = _2Pi / (Real) N;

        cudaMalloc(&cuda_W, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_Y, N*sizeof(cuFloatComplex));
        W = (Comp*) malloc(N*sizeof(Comp));

        init_w_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(cuda_W, w, N);
        cudaDeviceSynchronize();
        cudaMemcpy(W, cuda_W, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }

    ~DifDitFFT_cuda(){
        free(W);
        cudaFree(cuda_W);
        cudaFree(d_Y);
    }

    void dft(Comp* Y) {
        using namespace std;

        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        double tt = omp_get_wtime();
        for(int i = L - 1; i >= 0; --i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditdiffft_kernel<<<gridDim, blockDim>>>(d_Y, cuda_W, i, N);
        }

        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }

    void idft(Comp* Y) {
        using namespace std;

        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);

        double tt = omp_get_wtime();
        for(int i = 0; i < L; ++i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditdifidft_kernel<<<gridDim, blockDim>>>(d_Y, cuda_W, i, N);
        }
        complex_div_real<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_Y, N);
        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }
};


// cuda bluesteinFFT
struct BluesteinFFT_cuda {
    static constexpr char Name[] = "BluesteinFFT_cuda";
    const std::size_t N;
    const bool pow2;
    Comp* A; 
    Comp* B; 
    Comp* C;
    cuFloatComplex* cuda_C;
    cuFloatComplex* d_X;
    cuFloatComplex* d_Y;
    cuFloatComplex* d_A;
    cuFloatComplex* d_B;
    double tt, run_time = 0;

    union {
        char c;
        DitFFT_cuda dit;
        DifDitFFT_cuda difdit;
    };

    BluesteinFFT_cuda(std::size_t N) : N(N), pow2(isPow2(N)) {
        using namespace std;
        if (pow2) {
            new(&dit) DitFFT_cuda(N);
            return;
        }

        new(&difdit) DifDitFFT_cuda(ceil2(N + N - 1));
        auto w = Pi / (Real) N;
        // allocate memory in GPU
        cudaMalloc(&d_X, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_Y, N*sizeof(cuFloatComplex));
        cudaMalloc(&d_A, difdit.N*sizeof(cuFloatComplex));
        cudaMalloc(&d_B, difdit.N*sizeof(cuFloatComplex));
        cudaMalloc(&cuda_C, N*sizeof(cuFloatComplex));

        // allocate memory for computation
        A = (Comp*) malloc(difdit.N*sizeof(Comp));
        B = (Comp*) malloc(difdit.N*sizeof(Comp));
        C = (Comp*) malloc(N*sizeof(Comp));
        
        // initialize C (W)
        init_C_kernel<<<(N + BLOCKSIZE - 1)/ BLOCKSIZE, BLOCKSIZE>>>(cuda_C, w, N);
        cudaMemcpy(C, cuda_C, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);



        B[0] = C[0];

        for (size_t i = 1; i < N; ++i){
            B[i] = B[difdit.N - i] = C[i];
        }
        for(size_t i = N; i < difdit.N - N + 1; ++i){
            B[i] = Comp{};
        }
        difdit.dft(B);
        cudaMemcpy(d_B, B, difdit.N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
    }

    ~BluesteinFFT_cuda() {
        if (pow2) {
            dit.~DitFFT_cuda();
            return;
        }
        difdit.~DifDitFFT_cuda();
        cudaDeviceReset();
    }

    void dft(Comp* Y, const Comp* X) {


        if (pow2) {


            dit.dft(Y, X);
            return;
        }
        using namespace std;


        cudaMemcpy(d_X, X, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);

        
        tt = omp_get_wtime();
        complex_mul_conj_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_X, cuda_C, N, difdit.N);


        cudaMemcpy(C, cuda_C, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);// ☺☺☺☺☺☺☺☺☺☺☺



        // perform dft on A
        for(int i = difdit.L - 1; i >= 0; --i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((difdit.N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditdiffft_kernel<<<gridDim, blockDim>>>(d_A, difdit.cuda_W, i, difdit.N);
        }
        
        complex_self_mul_kernel<<<(difdit.N + BLOCKSIZE) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_B, difdit.N);
        // perform idft on A
        for(int i = 0; i < difdit.L; ++i){
            auto h = size_t{1} << i;
            dim3 blockDim(h < 1024 ? (BLOCKSIZE + h - 1) / h : 1, h < 1024 ? h : 1024);
            dim3 gridDim((difdit.N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.y - 1) / blockDim.y);
            ditdifidft_kernel<<<gridDim, blockDim>>>(d_A, difdit.cuda_W, i, difdit.N);
        }
        
        complex_div_real_mul_conj<<<(difdit.N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_Y, cuda_C, difdit.N, N);

        cudaDeviceSynchronize();
        
        run_time += omp_get_wtime() - tt;
        tt = omp_get_wtime();

        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }

    void idft(Comp* Y, const Comp* X) {
        if (pow2) {
            dit.idft(Y, X);
            return;
        }
        using namespace std;

        cudaMemcpy(d_X, X, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);
        cudaMemcpy(d_Y, Y, N*sizeof(cuFloatComplex), cudaMemcpyHostToDevice);


        complex_mul_kernel<<<(N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_X, cuda_C, N, difdit.N);
    
        dim3 blockDim(32, 32);
        for(int i = difdit.L - 1; i >= 0; --i){
            auto h = size_t{1} << i;
            dim3 gridDim((difdit.N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.x - 1) / blockDim.x);
            ditdiffft_kernel<<<gridDim, blockDim>>>(d_A, difdit.cuda_W, i, difdit.N);
        }

        complex_self_conj_mul_kernel<<<(difdit.N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_B, difdit.N);

        for(int i = 0; i < difdit.L; ++i){
            auto h = size_t{1} << i;
            dim3 gridDim((difdit.N/(2*h) + blockDim.x - 1) / blockDim.x, (h + blockDim.x - 1) / blockDim.x);
            ditdifidft_kernel<<<gridDim, blockDim>>>(d_A, difdit.cuda_W, i, difdit.N);
        }
        complex_div_real_mul<<<(difdit.N + BLOCKSIZE - 1) / BLOCKSIZE, BLOCKSIZE>>>(d_A, d_Y, cuda_C, difdit.N, N);
        cudaMemcpy(Y, d_Y, N*sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);
    }

    

};



void transComplex2DoubleArray(Comp* complex, Real data[],size_t N){

        for (int i = 0 ; i < N ; i++){
    
            data[2 * i] = complex[i].real();
            data[2 * i + 1] = complex[i].imag();
    
        }
    }
    
    void transDoubleArray2Complex(Comp* complex, Real data[], size_t N){
    

        for (int i = 0 ; i < N ; i++){
    
              complex[i].real( data[2 * i]);
              complex[i].imag( data[2 * i + 1]) ;
            
        }
    
    }
    

    

/*
    Function to run fft and report the time
    
    Params:
    ~ std::size_t sign: 1 represents forward dft, 0 represents inverse dft (idft)
    ~ std::size_t N: size of input data
    ~ Comp* input: input data
*/
void fft_cuda_1d(std::size_t sign, std::size_t N, Comp* input, Comp* output){
    using namespace std;
    cfft_plan plan = make_cfft_plan (N);
    Real data[N*2];
    // trans the Comp input date to data in pocketfft
    transComplex2DoubleArray(input,data,N);
    // neither the cfft_plan, use pocketfft to calculate the final result
     if( sign == 1){
            cfft_forward (plan, data, 1.);
        } else{
            cfft_backward (plan, data, 1./N);
     }
     // trans back the data in pocketfft to Comp output
    transDoubleArray2Complex(output,data,N);
    destroy_cfft_plan (plan);

}


/*
    Function to run fft and report the time
    this version2 let the make_plan run in advance and pass the plan as a param 
    Params:
    ~ std::size_t sign: 1 represents forward dft, 0 represents inverse dft (idft)
    ~ std::size_t N: size of input data
    ~ Comp* input: input data
*/
void fft_cuda_1d_v2(cfft_plan plan, std::size_t sign, std::size_t N, Comp* input, Comp* output){
    using namespace std;
    Real data[N*2];
    // trans the Comp input date to data in pocketfft
    transComplex2DoubleArray(input,data,N);
    // neither the cfft_plan, use pocketfft to calculate the final result
     if( sign == 1){
            cfft_forward (plan, data, 1.);
        } else{
            cfft_backward (plan, data, 1./N);
     }
     // trans back the data in pocketfft to Comp output
    transDoubleArray2Complex(output,data,N);

}
/**
    axis1,axis2 is (1,2) (1,3) (2,3) 
*/
template<class T,class T1>
void swapaxes(T* a, T1* b, int depth, int rows, int cols, int axes1, int axes2)
{
    for (int k = 0; k < depth; k++)
    {
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < cols; j++)
            {
                if (axes1 == 1 && axes2 == 2  || axes1 == 2 && axes2 == 1)
                    b[i * depth * cols + k * cols + j] = a[k * rows * cols + i * cols + j];
                else if (axes1 == 1 && axes2 == 3  || axes1 == 3 && axes2 == 1)
                    b[j * rows * depth + i * depth + k] = a[k * rows * cols + i * cols + j];
                else if (axes1 == 2 && axes2 == 3  || axes1 == 3 && axes2 == 2)
                    b[k * rows * cols + j * rows + i] = a[k * rows * cols + i * cols + j];

            }
        }
    }
    return;
}





template<class T,class T1>
void transpose(T* a, T1* b,int rows,int cols)
{
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
            b[j * rows + i] = a[i * cols + j];
		}
	}
	return;
}


/*
    Function to run fft2d
    
    Params:
    ~ std::size_t sign: 1 represents forward dft, 0 represents inverse dft (idft)
    ~ std::size_t N: size of input data
    ~ Comp* input: input data
*/
void fft_cuda_2d(std::size_t sign, std::size_t H ,std::size_t W , Comp* input, Comp* output){
    using namespace std;
    Comp* temp = (Comp*) malloc( H * W * sizeof(Comp) );
    for (int i = 0 ; i < H ; i++) {
        fft_cuda_1d(sign, W ,input + i * W , temp + i * W);
    }
    // transpose the matirx
    Comp* temp_transpose = (Comp*) malloc( H * W * sizeof(Comp) );
    transpose(temp,temp_transpose,H,W);
    Comp* temp_output = (Comp*) malloc( H * W * sizeof(Comp) );
    for (int i = 0 ; i < W ; i++) {
        fft_cuda_1d(sign, H ,temp_transpose + i * H , temp_output + i * H);
    }
    // transpose the matirx
    transpose(temp_output,output,W,H);
    free(temp);
    free(temp_transpose);
    free(temp_output);
}

/*
    Function to run fft2d
    invoke fft_cuda_1d_v1() to subtitute origin fft_cuda_1d()    
    Params:
    ~ std::size_t sign: 1 represents forward dft, 0 represents inverse dft (idft)
    ~ std::size_t N: size of input data
    ~ Comp* input: input data
*/
void fft_cuda_2d_v2(cfft_plan plan_h, cfft_plan plan_w,  std::size_t sign, std::size_t H ,std::size_t W , Comp* input, Comp* output){
    using namespace std;
    Comp* temp = (Comp*) malloc( H * W * sizeof(Comp) );

    clock_t startTime,endTime;
    startTime = clock(); 
    for (int i = 0 ; i < H ; i++) {
        fft_cuda_1d_v2(plan_w, sign, W ,input + i * W , temp + i * W);
    }
    endTime = clock();
    // transpose the matirx
    Comp* temp_transpose = (Comp*) malloc( H * W * sizeof(Comp) );
    transpose(temp,temp_transpose,H,W);
    Comp* temp_output = (Comp*) malloc( H * W * sizeof(Comp) );

    startTime = clock(); 
    for (int i = 0 ; i < W ; i++) {
        fft_cuda_1d_v2(plan_h, sign, H ,temp_transpose + i * H , temp_output + i * H);
    }
    endTime = clock();
    // transpose the matirx
    transpose(temp_output,output,W,H);
    free(temp);
    free(temp_transpose);
    free(temp_output);
}
/*
    Function to run fft3d
    

    copy from _raw_fft in numpy

    Params:
    ~ std::size_t sign: 1 represents forward dft, 0 represents inverse dft (idft)
    ~ std::size_t N: size of input data
    ~ Comp* input: input data
*/
void fft_cuda_3d(std::size_t sign, std::size_t D , std::size_t H ,std::size_t W , Comp* input, Comp* output){
    // original code of fft3d ,aim to implement a general fft function for 2d and 3d ,even for 1d.

    using namespace std;
    
    cout << "the first fft_cuda_1d" << endl;
    
    for (int i = 0 ; i < D * H ; i++){
        fft_cuda_1d(sign, W, input + i * W, output + i * W );
	for (int j = 0 ; j < W ; j++){
	//   cout << output[ i * W + j ].real() << "," ;
	}
	// cout << endl;
    }
    

    cout << "the second fft_cuda_1d" << endl;
    swapaxes(output,input,D,H,W,2,3);
    for (int i = 0 ; i < D * W ; i++){
        fft_cuda_1d(sign, H, input + i * H, output + i * H );
	for (int j = 0 ; j < W ; j++){
	  // cout << output[ i * W + j ].real() << "," ;
	}
	// cout << endl;
    }
    swapaxes(output,input,D,W,H,2,3);

    cout << "the third fft_cuda_1d" << endl;
    swapaxes(input,output,D,H,W,1,3);
    for (int i = 0 ; i < W * H ; i++){
        fft_cuda_1d(sign, D, output + i * D, input + i * D );
	for (int j = 0 ; j < W ; j++){
	  //  cout << input[ i * W + j ].real() << "," ;
	}
	// cout << endl;
    }
    swapaxes(input,output,W,H,D,1,3);

}

void fft_cuda_3d_v2(cfft_plan plan_h, cfft_plan plan_w, cfft_plan plan_d, std::size_t sign, std::size_t D , std::size_t H ,std::size_t W , Comp* input, Comp* output){
    // original code of fft3d ,aim to implement a general fft function for 2d and 3d ,even for 1d.

    using namespace std;
    
    for (int i = 0 ; i < D * H ; i++){
        fft_cuda_1d_v2(plan_w, sign, W, input + i * W, output + i * W );
    }

    swapaxes(output,input,D,H,W,2,3);
    for (int i = 0 ; i < D * W ; i++){
        fft_cuda_1d_v2(plan_h, sign, H, input + i * H, output + i * H );
    }
    swapaxes(output,input,D,W,H,2,3);

    swapaxes(input,output,D,H,W,1,3);
    for (int i = 0 ; i < W * H ; i++){
        fft_cuda_1d_v2(plan_d, sign, D, output + i * D, input + i * D );
    }
    swapaxes(input,output,W,H,D,1,3);
}
/*
    Function to run benchmark (cufft)
    
    Params:
    ~ Comp* Data: input data
    ~ int LENGTH: length of input data
*/
void benchmark(Comp* Data, int LENGTH){
    cufftDoubleComplex* CompData = (cufftDoubleComplex*) malloc(LENGTH*sizeof(cufftDoubleComplex));
    for(int i = 0; i < LENGTH; ++i){
        CompData[i].x = real(Data[i]);
        CompData[i].y = imag(Data[i]);
    }

    cufftDoubleComplex *d_fftData;
    cudaMalloc((void**)&d_fftData, LENGTH * sizeof(cufftDoubleComplex));// allocate memory for the data in device
    cudaMemcpy(d_fftData, CompData, LENGTH * sizeof(cufftDoubleComplex), cudaMemcpyHostToDevice);// copy data from host to device

    cufftHandle plan;// cuda library function handle
    cufftPlan1d(&plan, LENGTH, CUFFT_Z2Z, 1);//declaration
    // cufftPlan1d(&plan, LENGTH, CUFFT_C2C, 1);//declaration
    
    double tt = omp_get_wtime();
    cufftExecZ2Z(plan, (cufftDoubleComplex*)d_fftData, (cufftDoubleComplex*)d_fftData, CUFFT_FORWARD);//execute
    cudaDeviceSynchronize();//wait to be done

    cudaMemcpy(CompData, d_fftData, LENGTH * sizeof(cufftDoubleComplex), cudaMemcpyDeviceToHost);// copy the result from device to host

    Check_CUDA_Error("error");

    cufftDestroy(plan);
    cudaFree(d_fftData);
    for(int i = 0; i < LENGTH; ++i)
        Data[i] = Comp(CompData[i].x, CompData[i].y);

    free(CompData);
    
}


void cufft_1d_float(Comp* Data, int LENGTH, std::size_t sign){
    cuFloatComplex* CompData = (cuFloatComplex*) malloc(LENGTH*sizeof(cuFloatComplex));
        for(int i = 0; i < LENGTH; ++i){
            CompData[i].x = real(Data[i]);
            CompData[i].y = imag(Data[i]);
        }

    cuFloatComplex *d_fftData;
    cudaMalloc((void**)&d_fftData, LENGTH * sizeof(cuFloatComplex));// allocate memory for the data in device
    cudaMemcpy(d_fftData, CompData, LENGTH * sizeof(cuFloatComplex), cudaMemcpyHostToDevice);// copy data from host to device

    cufftHandle plan,plan_inverse;// cuda library function handle
    cufftPlan1d(&plan, LENGTH, CUFFT_C2C, 1);//declaration
    cufftPlan1d(&plan_inverse, LENGTH, CUFFT_C2C, 1);//declaration
    
    double tt = omp_get_wtime();


    if (sign == 1){
        cufftExecC2C(plan, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_FORWARD);//execute
    }else {
        cufftExecC2C(plan_inverse, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_INVERSE);//execute
    }
    
    // cudaDeviceSynchronize();//wait to be done

    cudaMemcpy(CompData, d_fftData, LENGTH * sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);// copy the result from device to host

    Check_CUDA_Error("error");

    cufftDestroy(plan);
    cufftDestroy(plan_inverse);
    cudaFree(d_fftData);
    for(int i = 0; i < LENGTH; ++i)
        if(sign == 1)
            Data[i] = Comp(CompData[i].x, CompData[i].y);
        else
            Data[i] = Comp(CompData[i].x/LENGTH, CompData[i].y/LENGTH);

    free(CompData);

}



void cufft_2d_float(Comp* Data, int H, int W, std::size_t sign){
    cuFloatComplex* CompData = (cuFloatComplex*) malloc(H * W*sizeof(cuFloatComplex));
        for(int i = 0; i < H * W; ++i){
            CompData[i].x = real(Data[i]);
            CompData[i].y = imag(Data[i]);
        }

    cuFloatComplex *d_fftData;
    cudaMalloc((void**)&d_fftData, H * W * sizeof(cuFloatComplex));// allocate memory for the data in device
    cudaMemcpy(d_fftData, CompData, H * W * sizeof(cuFloatComplex), cudaMemcpyHostToDevice);// copy data from host to device

    cufftHandle plan,plan_inverse;// cuda library function handle
    cufftPlan2d(&plan, H , W, CUFFT_C2C);//declaration
    cufftPlan2d(&plan_inverse, H , W, CUFFT_C2C);//declaration
    
    double tt = omp_get_wtime();


    if (sign == 1){
        cufftExecC2C(plan, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_FORWARD);//execute
    }else {
        cufftExecC2C(plan_inverse, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_INVERSE);//execute
    }
    
    // cudaDeviceSynchronize();//wait to be done

    cudaMemcpy(CompData, d_fftData, H * W * sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);// copy the result from device to host

    Check_CUDA_Error("error");

    cufftDestroy(plan);
    cufftDestroy(plan_inverse);
    cudaFree(d_fftData);
    for(int i = 0; i < H * W; ++i)
        if(sign == 1)
            Data[i] = Comp(CompData[i].x, CompData[i].y);
        else
            Data[i] = Comp(CompData[i].x/H * W, CompData[i].y/H * W);

    free(CompData);

}

void cufft_3d_float(Comp* Data, int D ,int H, int W, std::size_t sign){
    cuFloatComplex* CompData = (cuFloatComplex*) malloc(D * H * W*sizeof(cuFloatComplex));
        for(int i = 0; i < D * H * W; ++i){
            CompData[i].x = real(Data[i]);
            CompData[i].y = imag(Data[i]);
        }

    cuFloatComplex *d_fftData;
    cudaMalloc((void**)&d_fftData, D * H * W * sizeof(cuFloatComplex));// allocate memory for the data in device
    cudaMemcpy(d_fftData, CompData, D * H * W * sizeof(cuFloatComplex), cudaMemcpyHostToDevice);// copy data from host to device

    cufftHandle plan,plan_inverse;// cuda library function handle
    
    cufftPlan3d(&plan, D, H, W, CUFFT_C2C);//declaration
    cufftPlan3d(&plan_inverse, D, H, W, CUFFT_C2C);//declaration
    std::cout << "cufftPlan3d end" << std::endl;
    double tt = omp_get_wtime();


    if (sign == 1){
        cufftExecC2C(plan, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_FORWARD);//execute
    }else {
        cufftExecC2C(plan_inverse, (cuFloatComplex*)d_fftData, (cuFloatComplex*)d_fftData, CUFFT_INVERSE);//execute
    }
    std::cout << "cufftExecC2C end" << std::endl;
    // cudaDeviceSynchronize();//wait to be done

    cudaMemcpy(CompData, d_fftData, D * H * W * sizeof(cuFloatComplex), cudaMemcpyDeviceToHost);// copy the result from device to host

    Check_CUDA_Error("error");

    std::cout << "cudaMemcpy end" << std::endl;
    cufftDestroy(plan);
    cufftDestroy(plan_inverse);
    cudaFree(d_fftData);
    for(int i = 0; i < D * H * W; ++i)
        if(sign == 1)
            Data[i] = Comp(CompData[i].x, CompData[i].y);
        else
            Data[i] = Comp(CompData[i].x/D * H * W, CompData[i].y/D * H * W);

    std::cout << "cudaFree end" << std::endl;

    free(CompData);

}


